FILE_LOCATION   =	"/ifs/scratch/c2b2/ip_lab/ddd2124"


setwd(FILE_LOCATION)


PC_FILE<-read.table("parchild.lst",header=F,stringsAsFactors=F)


poplist<-read.table("PopIndex.txt",header=F, stringsAsFactors=F)
poplist<-poplist[,1]


#enact these once you have the lengths and numbers files
Lengths<-read.table("Lengths.txt",header=F)
Numbers<-read.table("Numbers.txt",header=F)

pcl<-integer(0)
pcn<-integer(0)

for (i in seq(length(PC_FILE[,1])))
{
	if ((A<-which(poplist==PC_FILE[i,1]))<(B<-which(poplist==PC_FILE[i,2])))
	{pcl<-c(pcl,Lengths[A,B])
	 pcn<-c(pcn,Numbers[A,B])}
	else 	{pcl<-c(pcl,Lengths[B,A])
		 pcn<-c(pcn,Numbers[B,A])}
}

pclcdf<-ecdf(pcl)
pcncdf<-ecdf(pcn)

g<-length(Lengths[,1])



#####################################################
##################PC MATRICES########################

PCNMAT <- Lengths
for (i in seq(g-1))
{
	PCNMAT[i,] <- pcncdf(Numbers[i,])
	PCNMAT[i,] <- apply(rbind(PCNMAT[i,],1-PCNMAT[i,]), 2, min)
}

PCLMAT <- PCNMAT
for (i in seq(g-1))
{
	PCLMAT[i,] <- pclcdf(Lengths[i,])
	PCLMAT[i,] <- apply(rbind(PCLMAT[i,],1-PCLMAT[i,]), 2, min)
}


PCNMAT2 <- PCNMAT + t(PCNMAT)
write.table(PCNMAT2,"PCNMAT.txt",row.names=F,col.names=F,quote=F)
rm(PCNMAT2, PCNMAT)
gc()
PCLMAT2 <- PCLMAT+ t(PCLMAT)
write.table(PCLMAT2,"PCLMAT.txt",row.names=F,col.names=F,quote=F)
rm(PCLMAT2, PCLMAT)
gc()

###################average length

AverageLengths <- Numbers / (Lengths + .001)
pca <-integer(0)

g<-length(AverageLengths[,1])

for (i in seq(length(PC_FILE[,1])))
{
	if ((A<-which(poplist==PC_FILE[i,1]))<(B<-which(poplist==PC_FILE[i,2])))
	{pca<-c(pca,AverageLengths[A,B])}
	else 	{pca<-c(pca,AverageLengths[B,A])}
}

pcacdf<-ecdf(pca)

PCAMAT <- AverageLengths
for (i in seq(g-1))
{
	PCAMAT[i,] <- pcacdf(AverageLengths[i,])
	PCAMAT[i,] <- apply(rbind(PCAMAT[i,],1-PCAMAT[i,]), 2, min)
}

PCAMAT <- PCAMAT + t(PCAMAT)
write.table(PCAMAT,"PCAMAT.txt",row.names=F,col.names=F,quote=F)
rm(PCAMAT, pcacdf, pca)
gc()
################

quit()


